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ABSTRACT 

We present a new equation of state for infinite systems (symmetric, asymmetric and neutron matter) based on an extended Skyrme 
functional constrained by microscopic Brueckner-Bethe-Goldstone results. The resulting equation of state reproduces with very good 
accuracy the main features of microscopic calculations and it is compatible with recent measurements of two times Solar-mass neutron 
stars. We provide all necessary analytical expressions to facilitate a quick numerical implementation of quantities of astrophysical 
interest. 
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1. Introduction 


A key ingredient for many astrophysical calculations is a reliable Equation of State (EoS) for isospin asymmetric matter, covering 
from symmetric nuclear matter (SNM) to pure neutron matter (NM), from low to high densities (^4-5 times saturation density). 
In this respect, a popular application is by instance the description of neutron stars (NS) properties as the mass-radius relation or 
the ir inhomogeneous c rust. Restricting ourselves to the category of nucleonic EoS, one of the most popular EoS is the one derived 
by iBaldo et^ (Il997l) . It has been obtained within the context of Brueckner-Bethe-Goldstone (BBG) many-body theory using the 
Argonne vl4 potential plus the Urbana model for the three-body nuclear interaction. Such an EoS has been tabulated for given 
values of the density of the system. Eor su ch a reason , it is c ustomary to fit the EoS with some analytical expressions which are 
much simple to handle in numerical codes (iTvrrel et alJ (l2013h ) suited for astroph ysical simulatio ns. 

A possible alternative to the fit is the use of an effective Skyrme interaction (ISkvrmd (Il959l) ) as early suggested by ICao et all 
i2006l) . These authors have shown that it is possible to fit an effective Skyrme functional (that they named LNS) on BBG, conserving 
some of the main featu res of the original BBG EoS . The main advantage of using a functional instead of a generic interpolation, 
as done for instance bv iHaensel & PotekhinI (|2003), is that once the parameters of the functional are fixed, all basic properties as 
pression or symmetry energy can be simply obtained by standard derivative operations. In the case where the vector part of the 
functional is also taken into account, as for the case of a functional derived from a complete Skyrme interaction, with a simple 
formalism based on the Linear Response theory iPastore et al.l (l2015l) . one can also describe collective phenomena within the NS. 
The la tter play a crucial role in descri bing different phenomena in the NS as the thermal properties of the inner crust (IChamel et al.l 
i2013h ) or the neutrino mean free path llwamoto & PethickI (Il982h. Another importan t advan tage is that the same functional can also 
be consistently used to describe the region of the crust of the NS dChamel & Haensell (l2008h ). thus allowing for a unified description 
of the star. Although LNS gives a nice reproduction on infinite matter properties up to two times saturation density, its EoS in pure 
neutron matter (PNM) remarkably deviates from the BBG results, leading to a different behavior of the s ymmetry energy at high 
density. It follows that the LNS EoS supports only NS with mass lower than 1.6 Solar-masses as shown bv ISingh et al.l (l2013h . The 
authors of (iGambacurta et akl (1201 ih ) have recently refitted the LNS functional, but hte new LNSl and LNS5 do not substantially 
improve the properties of the homogeneous nuclear medium as compare d to the original L NS, although they improve the description 
of finite nuclei. In the present article, we generalize the analysis done bv ICao et al.l (l2006l) concerning the possibility of constraining 
a phenomenological S kyrme functional on microscopic results, by using an extended Skyrme functional which includes up to 6th 
order deriyatiye terms (ICarlsson et akl (l2008h : lRaimondi et al.l (1201 Ih ). aiming at giying a reliable EoS also in the high density region 
and thus in better agreement with BBG results. _ 

Indeed, the Skyrme interaction can be interpreted as a low-momentum expansio n of a finite-range interaction (ISkyrmel(ll959h ). 
The standard form of the interaction is the one giyen by IVautherin & BrinM (Il972h and it takes into account only gradients terms 
up to the second pow er, as for LNS. Although this can be yiewed as a good approximation to be used in finite nuclei calculations 
(iBender et al.1 (l2003h ). it is not adapted to the study of dense nuclear matter. For example, the standard Skyrme interaction is not 
able to reproduce at the same time the correct isoyector splitting of the effectiye mass of BBG results (iBaldo et al.l (l2014bl) ) and 
the high density behayior of nuclear matter. Among the different strategies one can adopt to oyercome these difficulties, the most 
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promising implies either the addition of extra de nsity dependencies on the velocity dependent terms (IChamel & Goriel^ (1201 (ih i 
or the inclusion of higher order derivative terms (ICarlsson et al.l (l2008h i. We prefer to follow the latter approach sin ce it allows to 
grasp t he correct behavior of the EoS of BBG calculations especially at high density. In particular, it has been shown (iDavesne et al.l 
(l2015h i that the different terms follow a precise hierarchy and thus high order terms have stronger influence on the high density 
pa rt, inducing almost negligible mo dification to the low density part. Such a result is also in good agreement with previous findings 
of ICarlsson & Dobaczewsld (l2010ll . based on Density Matrix Expansion methods in finite nuclei. To this respect, the functional 
presented here can be considered as the natural extension of the LNS one to correct the high density region. In the present article, 
we thus present an extended Skyrme functional, hereafter called LYVAl, for a proper treatment of these higher order gradients, 
giving all necessary analytical expressions for astrophysical calculations. A numerical code and the tabulated values of this new 
EoS will be available at the CompOSE webpageQ. 

The article is organized as follows. In Sec.|2]we present the general formalism of the extended Skyrme functional. In Sec.|3]we 
give the general formula for the binding energy per particle for isospin asymmetric nuclear matter, while in Sec. |4]we present the 
case of polarized matter. In Sec.|5] we study the behavior of the symmetry energy. In Secl6]we further discuss the behavior of the 
effective mass and in Sec.|2]we examine the applications of our model to the description of a NS. Our conclusions are then given in 
Secg] 


2. Extended Skyrme interaction 

The m ost general form of the Skyrme functional up to 6th order in the gradient expansion has been derived by ICarlsson et al.l 
(l2008h . In the present a rticle, we prefe r to rel ate this functional to an effective interaction, thu s reducing the number o f free coupling 
constan ts, as shown in (IDavesne et al.l (l2013h i. The corresponding Skyrme interaction reads (iRaimondi et al.l (1201 ih i lPavesne et al.l 
(I2014blri 


= 4°' (l + (l + xfp^) [k'2 + k^] + tf (l + xfPa) k' • k + iff (l + xfp^) «“(R) 

+ iff (l + xf [(k'2 + k2)2 + 4(k' ■ k)2] + ff (l + xfp^) (k' ■ k)(k'2 + k^) 


-ff (l + xfp^) (k'2 + k^) [(k'2 + k2)2 + 12(k' ■ kf] + ff (l + xfPa) (k' ■ k) [3(k'2 + k^)^ + 4(k' ■ k)^] . 


( 1 ) 


The notations used here are standard and more details can be found in (iBender et al.l (l2003h l. The spin-orbit and tensor terms are 
here discard ed since they do not con tribute to the total EoS, although they do contribute to its multipolar partial wave decomposition, 
as shown bv IDavesne et al.l (l2014all . The corresponding functional form can be obtained by performing an average on Hartree-Eock 
states. Results for an homogeneous medium are given in Eq. (|5]l. An important advantage of deriving a functional from an effective 
interaction is that when applied to a single nucleon the functional leads to vanishing internal energy. This is not au t omati cally 
guarantee for a generic phenomenological functional thus implying that the nucleon can interacting with itself (IChamell (l2010h i. 
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et al. (l2015h . to sorne results of a BBG calculation dBaldo et al.l (Il997h : iBaldol (l2014lri . based on the microscopic Argonne vl4 
nucleon-nucleon two-body interaction pl us the Urbana mo del for the three-body term. Although more recent BBG calculations 
for the EoS are available (as for instance FZhou et akl (l2004l) i. the complete determinati on of the pa r amete rs requires also the SP- 
decomposition of the potential energy. To the best of our knowledge, the BBG results of (iBaldo et al.l (Il997h i are the most complete, 
also including results for ot her important a s trophy sical quantities. There exist other ah init io results obtained from chiral effectiv e 
field calculations at Viow-it (iHebeler et al.l (1201 ih l or the many-body perturbatio n theory ( Bogn er et al.l (l2005h : iRotivall (l2008h i. 
which could also be used to fix extended Skyrme parameters with a similar quality dPavesne et al.l ( 2015 11. However, they presently 
cover a na rrower densi t y rang e than BBG results, and are therefore not yet suited for our present purpose. We thus rely on the BBG 
results of (iBaldo et al.l (119970 1. The inputs for the fit include the projection of the energy per particle in the different spin (S ) and 
isospin (T) channels in symmetric nuclear matter (SNM), and the EoS of both SNM and pure neutron matter (PNM). As discussed 
bv IDavesne et al.l (12015h . no density-dependent term (i.e. tj^\ x®, a parameters) is required to get satisfactory fits, but the resulting 
parameterizations give too low a value for the Landau effective mass of SNM at saturation (m* jm ^ 0.4). In this paper, the density- 


dependent term is thus taken into account, and we have fixed its parameters to the values a = 1/6, tf^ - 13763 [MeVfm^"^®], 
and Xj*'^ = 0.3. In such a way, we can properly constrain the higher order derivative terms whose role is mainly to give the correct 
asymptotic behavior at high density. 

To fix the remaining parameters we have proceeded in two steps. In a first step, we considered the interaction ([TJ up to 4th order 
in the gradient expansion only. The parameters tf\ xf\ {n - 0,2,4) have been determined by fitting BBG results for SNM EoS and 
{S, T) channels. However, since the resulting EoS for PNM is too repulsive at high densities, we have added 6th order parameters 
on top of the previously determined 4th order values. To keep the quality of the SNM EoS we have imposed the values = 0, 
X® = 0, and x® = -5/4. The remaining 6th order parameter tf^ is then determined by fitting the PNM EoS in the full density 
interval. The resulting parameters of the LYVAl interaction are given in Table [T] 

In Eigure[Tl are displayed the EoS for SNM and PNM (panel a), and the SNM potential energy decomposition in the different 
spin-isospin {S, T) channels (panel b). Our fit is clearly very good for both EoS. The results obtained using LNS are also displayed 


* http://compose.obspm.fr/ 
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Table 1. Parameters of the extended LYVAl Skyrme interaction, with a = 1/6, = 13763 [MeVfm^’^“], and x® = 0.3. 
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Fig. 1. (Colors online) Equations of state of SNM and PNM (panel a) and projections (5, T) in SNM (panel b), both exp ressed in MeV. The solid 
lines represent the result obtained with our extended Skyrme interaction, while the dots represents the EoS obtained bv iBaldo et alj l ll997ll . The 
LNS results are represented by dashed-lines. 


(dashed lines), and one can see a rapid deviation of the PNM EoS, starting from n 0.4 fm“^; one can expect this deviation also 
manifests for other quantities as the symmetry ener gy at high values of t he density. Let us now turn to the results for the (S, T)- 
channels shown in panel b. As already discussed in ( Lesinski et al.l (l2006l) ). a general drawback of the standard Skyrme functional 
is that the simultaneous reproduction of the {S, 7’)-channels is very difficult, to say the least. In the figure one can see the particular 
LNS case, which fails to reproduce BBG results. In contrast, with the extended functional the S - I channels are nicely fitted in 
the full range of density values, whereas the S - 0 channels show a deviation for n > 0.6 fm“^ as a consequence of our giving 
more weigh to PNM data in the fit. All the other quantities presented hereafter in the article have not been fitted, and they can be 
considered as a prediction of our model. 


It is worth mentioning that there are some other functionals whi ch have been developed wit h particular attention to the properties 
of NS. Among the no n-relativistic ones, we consider the BCPM (iBaldo et al.l (l2013Ll^l4al) : ISharma et aT] (l2015h ) and the BSk 
family dGorielv et alJ (l2009l l2()T^ ). The BCPM functional has been derived in a complete Khon-Sham scheme, thus not related 
to any interaction, and it has been explicitly constrained to reproduce BBG results in homogeneous matter. The BSk models have 
been derived from an effective Skyrme interaction with the addition of a power of the density into the momentum dependent terms 
of the standard Skyrme interaction. The BSk model have been constrained on several nuclear observables as masses and radii of 
finite nuclei together with additional pseudo-observables of homogeneous nuclear matter. In Figl2l we compare the EoS in both 
SNM and PNM obtained with LYVAl, the BCPM functional and three representative BSk interactions, namely BSkl9, BSk20 and 
BSk21 (IChamel et alJ (1201 ih ). Since the BCPM functional fits the same microscopic EoS as LYVAl, we observe that the results are 
almost on top of each other, except in the saturation region where the BCPM has been adjusted to give a value of saturation density 
no - 0.16fm“^ and E/A = -16 MeV. We remind that the BCMP has been fitted up to n = 0.6 fm“^ and that beyond that value 
the microscopic results of iBurgio & Schulz^ (1201 Oh have been used. In order to be consistent, we will thus omit the points beyond 
0.6 fm“^ in this paper. 


Concerning the BSk functionals, we used the generalize d expressions given in (iLesinski et all (l2007h ) to obtain their S T decom 


the coupling constants in front of the so-called 
) to avoid the appearance of spurious ferromag- 


position. It is important to notice that as for the SLy4 cas e (IChabanat et"^ (Il997h) 

term are switched to zero. This choice is justified injJChamel_&_Goriely|j 20^ ____ ^ 

netic phase-transitions in the homogeneous medium (iMargueron et al.l (l2002h ) md anomalous behavior of the entropy. In FigOJ 
we compare the results obtain ed with the BSk mod els and the BBG calculations. We observe that the BSk behave better than any 
standard Skyrme interactions (ILesinski et alJ (l2007l) ) since in the low-density region (a; no) the BSk give the correct sign and trend 
of the energy per partic le. On the same figure we also report the chiral effective field theory (j^^-EFT) calculations at low momentum 
k (iHebeler et al.l (1201 Ih ). The (j^^-EFT) results are in very good agreement with the BBG results apart from the (S=1,T=1) channel. 
Such a comparison gives us the level of uncertainty related to the adopted interaction and/or calculation technique. 
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Fig. 2. (Colors online) Equations of state of SNM (panel a) and PNM (panel b) obtained using the LYVAl interaction (solid) the BSk models 
(dashed) and the BCMP functional (dotted). 



Fig. 3. (Colors online) Same as Fig[T] but for the BSk model discussed in the text. 


3. Energy per particle 

In this section we give the analytical expression of the binding energy per particle (E/A) for infinite systems with isospin unbalance, 
here called Asymmetric Nuclear Matter (ANM). Other relevant quantities can be easily derived from it, including SNM and PNM. 
It is convenient to define an isospin asymmetry parameter as 


Y = 


= 1 - 2Y„ 


( 2 ) 


where nn(p) is the neutron (proton) density, n - + n-g is the total density of the system and Yp is the proton fraction. When the 

asymmetry parameter is equal to F = 0 we are in the SNM case, while for F = 1 we are in the other extreme case, i.e. PNM. To 

present the expressions in a compact, yet transparent, form we define the coefficients a 
following functions of the asymmetry parameter 


(37rV2) 


1/3 


and b 


^ , and the 


F. - i[(l + F)-‘ + (l-F)T , 




(1-F)^] . 


The EoS in ANM reads 

+ i [3 _ (^2xf + 1)f 2] ^ [3 - (2x3 + 


EjA = 


5 2m 8 


80 

1 




1120 


K{5f- 


7/3 


+ 7F: 


) + Cf^(5FG7/3+7G2/3))]n^/3 


240 
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The constants are the following combinations of Skyrme parameters 

+ (5 + )4"\ 

cf - +1)4"^ + ( 24 "^ + i)ff. 


(6) 

(7) 


E/A [MeV] 



0 0.1 0.2 0.3 O.4., 0.5 0.6 0.7 0.8 

n 


Fig. 4. (Colors online) Equations of states in asymmetric nuclear matter as a function of the density n and asymmetry parameter Y. 


In Fig. m we show the binding energy per particle obtained with our new functional as a function of the density and the 



method to another. The goal of the present article is to prove that a simple Skyrme functional can grasp the main features of a more 
complicated calculation based on realistic nucleon-nucleon interaction, as a consequence we prefer not to do any fine tuning around 
the saturation point since it will not change the main conclusions of the present work. 

From Eq. |5j we can also extract other quantities, as the pressure of the system P - or the nuclear incompressibility 

K - 9n^ “ 9 If j as functions of the asymmetry parameter Y. In Tab|3 we report several relevant SNM quantities 

calculated at saturation density. Our parametrization gives a valu e of the incompres sibility of /T = 231 MeV at saturation density, 
which is within the range of acceptable values as discussed by iDutra et "^0 (l2012h . The third derivative of the EoS gives us the 
skewness Q. 


Table 2. Basic SNM properties calculated with the LYVAl parametrization given in Tab[T] the BCPM and the BSkl9-21 functionals at saturation 
density no- 



LYVAl 

BCPM 

BSkl9 

BSk20 

BSk21 

LNS 

no[fm^^] 

0.169 

0.160 

0.160 

0.160 

0.158 

0.175 

£/A[MeV] 

-17.02 

-16.00 

-16.08 

-16.80 

-16.05 

-15.31 

K[MeW] 

231 

214 

237 

241 

246 

211 

m* jm 

0.707 

1 

0.80 

0.80 

0.80 

0.825 

e[MeV] 

-463 

-881 

-298 

-282 

-274 

-384 

7[MeV] 

33.8 

31.9 

30.0 

30.0 

30.0 

33.4 

L[MeV] 

64.5 

53.0 

31.9 

37.4 

46.6 

61.5 

[MeV] 

-75.6 

-98.1 

-191.4 

-136.5 

-37.2 

-127.7 

e.vv«[MeV] 

464 

877 

473 

550 

710 

303 
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The two limiting cases of symmetric nuclear matter and pure neutron matter can be immediately obtained from Eq. 


13 


E/A 


SNM 


E/N 


PNM 


— —H-H—H- 

5 2m 8 ° 16 80 ° 280 ° 15 “ 

+ ^(1 - xf)Cnn + ^(1 - X3)f3nr‘ 

[C® + Cf>] + :Eb^ [C + Cf>] n^/3 + ^^6 [^(6) ^ ^(6)| ^3 ^ 


(8) 


(9) 


In Fig.|5] we show the evolution of the pressure P as a function of the density in SNM for LYVAl, BSkl9-21 and BCPM. The 
two areas represe nt the constra ints obtained on the EoS bv iDanielewicz et al.l (l2002h using experimental observations of heavy-ion 
collisions and bv lEuchsI (l2006h of experiments on kaons. We observe that the LYVAl functional is perfectly consistent with such 
results. 



Fig. 5. (Colors online) The Pressure as a function of the density in SNM for the different models considered in the article. See text for details. 

In Figl^ we compare the resulting EoS in PNM for the LYVAl functional and t he results of several microscopic calculations 
(lAkmal et^.l (ll998h : lLi & Schul^ (l2008h : lBaldo et af1(ll997h : lGandolfi et al.l (l2012h ). We observe that our EoS is compatible with 
these calculations up to three times the saturation density; beyond that value the different calculations strongly differ from each 
other. It is thus very important to adopt the same microscopic calculation to constrain both the PNM ans SNM EoS, otherwise, one 
would obtain non trustable results concerning the behavior of the symmetry energy. 



n [fm 

Fig- (Colors online ) The EoS in PNM for the different models discussed in the text. The shaded area represents the constraints extracted from 
dGandolfi etalld2012M . 


In Figlll we sho w the pressure P in PNM for the different models discussed here and some additional constraints derived in 
(iDanielewiczl (l2003h ) and based on the analysis of heavy ion collisions. We observe that both the BCPM and LYVAl models are in 
good agreement with the constraints extracted assuming a soft EoS. See (IDanielewiczl (l2003h ) for more details. 
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Fig. 7. (Colors online) The pressure as a functi on of the density in PNM for the different models considered in the article. The grey areas are 
derived from the analysis of dDanielewicd l l2003ft ). 


The sound ve l ocity v , of a system is obtained in the non-relativistic limit from the isothermal incompressibility dnIdP. We refer 
to (iHaensel et al.l (l20()7h ) for a more detailed discussion. The explicit expression for PNM reads 


m„c 


1 


-b^ [c® + C®] + -b^ + cf>] nl'^ + -b^ [c® + Cf ] 


( 10 ) 


In Fig|8] we show the evolution of Vj in PNM as a function of the density. We observe that our parametrization respects the 
causality principle (there is actually a maximum around 1.1 fm“^ with vjc - 0.97), and is thus reliable for the description of high 
density neutron matter. On the same figure we also compare the resu lts obtained w i th the BSk models. We notice that the BSk21 
model violates causality in PNM at n ^ 0.8 fm“^. As discussed in (iGorielv et al.l (1201 Oh ), the BSk models assure the causality 
principle inyS-equilibrium nuclear matter for the density ranges found in NS. 

From this analysis we exclude the BCPM model since its analytical expressions are strictly valid, by construction, only in the 
low density regime. 

We notice that although the causality principle is always respected by the LYVAl functional a speed of sound very close to the 
speed of light is clearly a symptom of the use of the non-relativistic approximation to treat matter very close to the relativisitic limit. 



Fig. 8. (Colors online) The speed of sound in PNM as a function of the density of the system. The notation is as in Fig[T] 


4. Polarized Matter 

A major drawback of several effective interactions is the presence of a spontaneous ferromagnetic transition (IVidaurre et al.l (Il984h ) 
at densities relevant for the physics of nuclei and or NS. However, such a spontaneous phase transition has not been observed so 
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far by any microscopic calculation (iPandharinande et^ (ll972h : lFantoni et^(l200lh : IVidana et al.l(l2002h '). It is thus interesting to 
determine the behavior of our interaction concerning this aspect. 

The expressions for the energy per particle in fully polarized pure neutron matter (PolPNM) reads 


E/A 


PolPNM 








+ (6)U6),^3 

15 ^ ^ 


( 11 ) 


where we have used the notatio n c — (67 tW^^. I n the following we compare the results obtained with the LYVAl functional and 
available BBG calculations of (iBombaci et al.l (l2006h l. More precisely, we show in Fig0 the difference of energy per particle 
between PolPNM and PNM, that is AE/A — E/A\„ - E’/AL,,., obtained with the LYVAl functional and the BBG results 

of (IBombaci et all (1200611 1. In order to be consistent, we consider only, for this particula r case, the re s ults up to ssS no since the 
treatment of the three-body term at high density is not the same used in BBG results of (iBaldo et al.l (Il997h l and used here for 
the ht of the LYVAl functional. We observe that the LYVAl as well as BSk20 follow pretty closely the BBG results, while the 
BSkl9(21) tends to underestimate (overestimate) the energy difference between the two systems. The LNS functional is not stable 
against polarization and at densities n 0.6 [fm“^] favors the appearance of polarized neutron matter. The BCPM functional has 
not been included in such analysis since the functional has not been tailored to describe polarized systems. 



Fig. 9. (Colors online) Energy difference between PolPNM and PNM for the different models considered in the text. 


5. Symmetry energies 

We give now the LYVAl expression for the isospin symmetry energy srin), which plays a crucial role in determining the composi¬ 
tion of the NS since the y6-equilibrium condition strongly depends on it. It follows that reproducing the symmetry energy not only 
at saturation, but also as a function of the density is a necessary condition to have a reliable extrapolation of the high density part 
of the NS. Starting from the complete expression of Eq. |5j we can expand the binding energy per particle up to second order in the 
following way 

j{n,Y)= j{n,0) + STin)Y^ + ... ( 12 ) 


to get the result 


Srin) — — a^n^^^ 
6m 




—a^ { rC® + C' 


< 2 ) 1 „ 5/3 


-(-— 
32 


-,(4) 


3 ' 


(4) I „7/3 


1 


48 




-( 6 ) 


(13) 


In Fig. [To] we show the evolution of the symmetry energy ej as a function of the density of the system. At saturation density, 
we obtain a value of the symmetry ene rgy srinp) — J - 3 3.8 MeV, a value compatible with most recent constraints on J obtained 
combining different experimental data (iTsang et al.l(l2009h ). Furthermore, we observe an excellent agreement up to several times the 
saturation density value between our results and the BBG ones. In the same figure, we also compare the evolution of the symmetry 
energy for the BCPM and BSk models. We observe that while BCPM gives by construction results which are essentially on top of 
ours for low density, the BSk give very different behaviors especially beyond saturation density. There is not agreement between 
diffetenMnictosco£ic_a 22 roaches_concerningJhe^ehavior_ofej;^e^on^^aturation_densit^jJVejeferjoJhe_discu^sionj^ 

et al. (I 2 OIOI) ). A possible way to figure out the corr ect trend of sj at high density is the predicted proton fraction and thus the 
possibility or not of allowing firect URCA process (iHaensell (Il995h ). We refer to SecjTjfor a more detailed discussion. We can 
anyhow anticipate that BSkl9-20 and LNS are not compatible with such additional constraints. In panel (b), we compare the results 
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at low density obtained with the different models and the constraints obtained in ( iDanielewicz & (l2014ll '). The large yellow 

area represents the constraints extracted by analyzing data on isobaric analog states (IAS), while the smaller area delimitated by the 
solid line contour has been obtained by studying data on the neutron skin properties of some selected nuclei. All the functionals 
considered in the text respect these constraints. 



Fig. 10. (Colors online) Symmetry energy as a function of the density of the system for different models (left panel). The dots correpond to the 
BBG calculations. In the right panel we compare the constraints extracted using IAS IDanielewicz & ( 120141) . See text for details. 


It is important to remind that a parabolic approximation was used in (iBaldo et al.l (11997)) to extract the symmetry en ergy and 
the validity of such an approximation has been tested only in the region n e fO. 0.41 fm~^ dBombaci & Lomb^dol (Il99lh: Gorielv 

et al. (l201f)l) ). In this case the symmetry energy is obtained as the difference between the EoS in PNM and SNM 





SNM 



PNM 


(14) 


In Figll II we compare the results for the LYVAl functional using the definition of Eq[T3]and EqlT4] We observe t hat usi ng the 
parabolic approximation we have a better reproduction of the high density part of the BBG results of (iBaldo et al.l (Il997l) ). The 
agreement is still not perfect due to the small overestimate of the EoS for PNM resulting in our fit and shown in Eig[T] We notice 
that in our calculations the two definitions of the symmetry energy of Eq[T3]and EqlT4]are essentially on top of each other up to 
11 ^ 0.6 [fm“^]. 



Fig. 11. (Colors online) Symmetry energy as a function of the density of the system using different definition of the symmetry energy. The dots 
represent the BBG results. See text for details. 


For completeness, we also define additional quantities related to the symmetry energy, as L — Ksym - and 

Qsym - 27»^ ^ . Their values at saturation density are reported in Tab.|2] The constraints on these quantities are much less strict, 

leading to larger error bars dPutra et al.l (l2012h ). However, the value obtained here by instance for L is compatible with some recent 
estimates (IChen et al.l (1201 ol) ) extracted from finite nuclei analysis. 


6. Effective mass 

The effective mass is directly related to some important processes as neutrino emissivity (lYakovlev et al.l (1200Ih ). We thus give the 
explicit expressions for the neutron (»?*) and proton (mp Landau effective masses at the Fermi surface as a function of the density 
and of the asymmetry of the system. For neutrons it reads 
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(15) 


The proton effective mass is simply obtained by replacing Y —> -Y in the previous expression. It is worth noticing that contrary to 
the standard Skyrme interaction, the effective ma ss for our pseudo-po tential has an explicit dependence on the momentum k, as it 
happens in the case of real BBG calculations. See lBecker et al.l (1201 4) for more details. 

In Fig. [ 12 ] we show the evolution of the effective mass for neutrons and protons as a function of the asymmetry parameter 
at saturation de nsity. We obs erve that the mass splitting Am - m* - m* has the correct sign and density behavior as compared 
to BBG results (iBaldo et alJ (l2014bll '). although the resulting effective masses are slightly lower at saturation in the SNM case 
im* jm ^ 0.7) compared to the BBG result {m*fm ^ 0.8). Such a difference can’t be further reduced by a better fine-tuning of the 
tf\ X® parameters without inducing side effects on other quantities. 
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Fig. 12. (Colors online) Neutron and proton e ffective mass at saturation density as a function of the asymmetry parameter Y. The dots correspond 
to the BBG calculations dBaldo et all ll2014blf ). the solid line to LYVAl interaction. 

In Fig.[T3](a), we compare the effective mass in SNM calculated with the LYVAl functional and the corresponding BBG results 
as a function of the density of the system. Although the difference between the two calculations increases with the density reaching 
at most 30% in the high density region, the asymptotic behavior is correct with a positive slope at high densities. On the same figure 
we also show the results obtained with the different BSk models (concerning BCPM, it has been fitted imposing the bare nucleon 
mass at all densities and all asymmetries). The major difference between the BSk and LYVA models is related to the high density 
behavior, where the former lead to a much smaller effective mass compared to BBG results. In Fig.[T3](b), we show the evolution of 
the difference between neutron and proton effective mass Am - m* /m - niplm at no = 0.169 fm“^ for the different models discussed 
in the present article. We notice that the BSkl9-20 models give a much larger splitting than the one observed with BBG calculations. 



Fig. 13. (Colors online) On the panel a, we show the neutron effective mass in SNM (Y=0) for the BBG calculations (dots), LYVAl (solid) and 
BSkl9-21 (dashed). On the panel h the evolution of the effective mass splitting at no = 0.169 fm“^ . 


The presence of magnetic fields inside a star could lead to spin-unbalanced systems. As discussed in SecHl such a configuration 
does not constitute the ground state of the system up to several times the saturation density, but the presence of an external field could 
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change the situation. Limiting ourselves to the case of polarized pure neutron matter, it is possible to write the explicit expression 
for the spin-up (t) and spin-down (X) effective mass as 

^ (. 6 ) 

+ \cfnlCp{l + + icf n^2^(l + 

o o 

+ [(1 + + (1 - A ) 5 / 3 ] + [(1 + - (1 - 

lo 10 

+ ^cf n4(l + A)4/3 + Icfnk^Ml + A)4/3 

-H n4(l + A)2^^[(l -H -H (1 - A)^/^] H- ?j 4(1 + A)2/3[(l + - (1 - 

+ ^Cfn4[{\ + + (1 - A)"^^] + -^Cf n4[(l + A)^'^ - (1 - A)"/^] 

lo lo 

where we have defined A = (p| - p[)lp. The coupling constants C^f’ are related to the if’, xf as follows 
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(18) 


The explicit expressions for the spin-down effective mass can be derived from EqlT6]by replacing A —> -A. In Fig. [T 4l w e show the 
evolution of the spin up (down) effective mass for the LYVAl model and the original BBG results (iBombaci et alJ ( 20061) 1 



Fig. 14. (Colors online) Spin up and spin down neut r on eff ective mass at saturation density as a function of the polarization parameter A. The dots 
correspond to the BBG calculations dBombaci et all l l2006h 'l. the solid line to LYVAl interaction. 

In FigfTSl we compare the effective mass for the spin up (down) component for the different functionals considered in the present 
article in polarized neutron matter. The BCPM results are not present here : this functional can not be used for polarized systems 
since the informations on the vector part are missing by construction. Contrary to the LYVAl functional, the BSk models do not 
produce any splitting in the effective mass : this is due to the absence of terms C^s • T (see Eqs. A1-A2 of iPastore et al.l (l2014tl ) in 
the funct ional which governs such a splitting. The effective mass given by BSkl9 is particularly high compared to BBG results. We 
refer to ( iGorielv et al.l (1201 Oh ) for a more detailed discussion. 

7. Neutron Stars 

In this section, we present the basic properties of a non-accreting NS at zero temperature using our interaction. To calculate the 
mass and the radius of a NS we have to solve the Tolman-Oppenheimer-Volkoff (TOV) equations for the total pressure P and the 
enclosed mass m 


dP{r) 

dr 

dm{r) 


Gm{r)e{r) 

yl 


P(r)g(r) \/ 47rr^P(r) y 
/ \ s{r)c^ / 


2Gm{r) 


dr 


47Tr^s(r) , 


(19) 
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Fig. 15. (Colors online) Same as Fi dl2l but for other models considered in the text. Up (down) triangle stands for spin up (down) component. 


where G is the gravitational constant and s{r) is the total energy density of the system. Since in our model we consider only neutrons, 
protons and electrons, the y6-equilibrium condition at each value of the density of the star translates into the equation fip+fi,, — Pe for 
the chemical potentials, the possible con tribution of mu ons being neglected. In Fig. [16] we show the proton fraction inside the star 
calculated using our Eo S. According to (Haense i (Il995h ). the direct URCA process, which is very important to have a fast cooling 
during stellar evolution (iLattimer et al.l (Il99lh ). can take place when the proton fraction is x ^ 0.11. We observe that our model 
predicts the possibility for direct URCA process at already 3 times saturation density. On the same plot, we also present the BBG 
results. The agreement between the two calculations is very good up to « ^ 0.6 fm“^. 

As anticipated in SecjSj the possibility of allowing or not a direct URCA process can be used to make some consideration 
concerning the behavior of the symmetry energy in the high density region. The yS-equilibrium condition for nucleonic equation of 
state can be related directly, within the parabolic approximation, to the symmetry energy, we have 


Pe = fin -fip~ 4 s 7'(«)(1 - 2Yp) ( 20 ) 

where fie,np, are the chemical potential of the different species included in the EoS. Erom EiglTb] we can conclude that only the 
LYVAl, BSk2i and BCPM functionals allow for a direct URCA process in NS. 



Fig. 16. (Colors online) Proton fraction Y,, as a function of the density of the star. The dashed line represents the proton fraction treshold to activate 
the direct URCA process during the cooling stage of the NS. The open dots represent the BBG calculations. 

To describe the structure of the NS we need to solve Eq. ([T9l l imposing the y6-equilibrium for each value of the density. In the 
most externalTaTCrejJhe_cru^tjUs^ossibleto.observeJhe^resenceofsttuctureSj,eithemiuclei,orjnoreexotic^aj^^ 

& Haensel ~(l2008h ). This part of the star will be described bv using the EoS of (iDouch in & Haensell (|200lh). which has been derived 
using the SLy4 functional and by means of the Compressible Liquid Drop Model JPouchin et al.l (l20()0l) ). It allows for a simple 
descriptionofbothJhe_cru^^innerand_outer}j^utaUoJor^heJi 3 uid;Cote^ansitioiuJ^ejnatch_ourEoS^ithJhe_OT 
& Haensel ~(l2001h ) at w 0.08 fm~^. This small inco nsistency in the EoS w ill not affect the value of the maximum mass of the star, 
but it introduces an error of at most 5% on its radius (iHebeler et al.l (l2010l) ). The study of the inhomogeneous phase of the NS with 
our functional will be the subject of a forthcoming study. Our results concerning the Mass-Radius relation are shown in Eig. ITTl a). 
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On the same plot, we also show the recent measurements of masses of NS (iDemorest et al.l (l2010h : lAntoniadis et alJ (l2013h '). which 
are both compatible for a 2Msun neutron star. We observe that our EoS gives a prediction compatible with the latest experimental 
measurements giving a maximal value of M = 1 . 96Ms,m in the non-rotating case. Th e inclusion of extra degrees of freedom as 
pions4^aonsorjl^£eronswouldaffec^ha^esultJ|Heiseibe^||^^iortOense^JSo0^i]J_as,^ell^,asJhe_effec^_ofrotationJSalgado 

et al. (ll994ah : IStergioulasl(l2003h l. In view of the results of (ISalgado et alJ(ll994b[) '). we could expect an increase of 10 - 15% for 
the results of our EoS. The detailed study of these effects goes beyond the scope of this paper and we leave it for the future. 



Fig. 17. (Colors online) In panel (a), we show the mass-radius relation for NS obtained with o ur EoS by solving TOV equations. The two horizontal 
bars refer s to the two recent of N S masses measurements = 1.97 ± 0.04 given in ( IDemorest et alj (120101) ) and = 2.01 ± 0.04 

given in dAntoniadis et al.l(l201.3h l. In panel (b), the mass-density relation. 

The radius of a NS is very difficult to extract from o b servations due to th e several hypothesis one has to do on the atmosphere 
composition. Different models (ISuleimanov et al.l (1201 ih i lSteiner et ^ (1201 Oh ) lead to slightly different values for the radius, but it 
is neve rtheless possible to give an upper value of around 12.5 km for a NS with a mass 1 AMsun- See also discussion in (lEortin et al.l 
(l2014h ). Erom the EoS of our functional, we get 11.6 km, in fair agreement with the original BBG results. In Eig. fTTl' b). we show the 
ev olution of the maxim um mass of the NS as a function of the central density of the star. It is worth noticing that recent constraints 
of (iKIahn et al.l (l2006h ) implies the absence of direct URCA process for NS within a mass range of 1 - 1 .SMsun- Erom Eigll7l(b). 
we can observe that the lowest value of the density at which URCA process take place correspond to a NS of mass 1.54 Mj„„. 


8. Conclusions 


We have presented a new nucleonic Equation of State based on the extended Skyrme functional. The inclusion of higher order 
derivative terms allow us to give a more precise description of the high-density region. By fixing the coupling constants of our 
functional on the EoS of ab initio calculations, we have shown the possibility of e xtracting analytically several quantities of strong 
astrophysical interest as incompressibility, pressure or effective mass. As shown in (iDavesne et al.l(l2015h ). the higher order gradients 
allow us to grasp the correct physical behavior obtained with a microscopic calculation. Our functional can be fitted to other 
microscopic calculations, thus providing a more powerful tool for astrophysics than a simple interpolation procedure. The Skyrme 
functional_^anbe_easiI^^m£lemente^to_£erforin£al£uIutionsjn_yHa^ere_offlie_stoj_jio£_onl^jnJfsjmifor^^ & 


Haensel ~(l2001 ): lPearson et al.l (l2012h ). 


WeJmvecorn£ared_ourresuIts^ithother_commonl^adoDtedJunctionalSjMmel^die^CPMjlBaldoeta^jl^^ 

et al. (l2015h ) and the BSk functionals ( iGorielv et al.l(l2009t) : IChamel & Gorielvl(l2010l) : lGorielv et al.l(l2010Ll2013h ). We have shown 
that our model is complementary to the results obtained to these different group since it aims at reproducing as accurately as possible 
BBG results in nuclear matter (including polarized matter), so that the functional can be used to describe both ground state properties 
and excited states. 

Due to the simplicity of the calculations, the formalism can be easily extended to properly include finite-temperature effects. In 
that case, it is necessa ry to replace the step func tion used to evalua te the Eo S of the s ystem b y a Eermi-Dirac distribution ( Bonche 
& Vautherin (Il98ll) ). These integrals can be approximated with analytical expressions as shown bv lAntial (Il993l) . The possible 
temperature depende nce of the coupling constants can be also studied by direct comparison with existing microscopic calculations 
at finite temperature dPandharinande & Ravenhalil (ll989l) : lL,eieune et"^ (Il986l) ). 
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